#/usr/bin/perl                                                                                                                           

# GLOBAL

my $grid_size_x = 300;
my $grid_size_y = 300;
my $locitot = 99;
my $name = "sim";

# Seed the random number generator
srand (0);

sub draw_rand{
    my ($n,$tot) = @_;
    my %h, @arr;
    while(keys(%h)<$n){
        $r = int(rand($tot));
        if(!$h{$r}){
            $h{$r}=1;
            push(@arr, $r);
        }
    }
    @arr = sort {$a <=> $b} @arr;
    return @arr;
}

# draw_rand_hash(n,tot) draws n random numbers from 0..(tot-1) without duplicates 
sub draw_rand_hash{
    my ($n,$tot) = @_;
    my %h;
    while(keys(%h)<$n){
        $r = int(rand($tot));
        if(!$h{$r}){
            $h{$r}=1;
        }
    }
    return %h;
}

# Draw random x cells, but since we have hex grid
# x = 0, 0.5, 1, 1.5 ... 511, 511.5
# No two cells have the same coordinates
sub draw_rand_cells{
    my ($n,$tot) = @_;
    my %h, $size=0;
    while($size<$n){
        $x = int(rand($tot));
	# either add 0.5 to both, or not
	$offset = 0.5*int(rand(2));
	$x+=$offset;
	if($offset<0.2){
	    $y = int(rand($tot/2))*1.5;	    
	} else {
	    $y = int(rand($tot/2))*1.5+0.75;	    
	}
        if(!$h{$x}{$y}){
            $h{$x}{$y}=1;
	    $size++;
        }
    }
    return %h;
}

# Start from sweep 1
# do not use my, so that these hashes are global
 %grid;
 %patterns;
 %indexx, %indexy; 
for $sw ($ARGV[0]..$ARGV[1]){

    my %grid_coords_x =();
    my %grid_coords_y =();
    open (in,"$name$sw/$name$sw.crypts");
    while(<in>){
	chomp;
	@row=split(/\,/);
	$grid_coords_x{$row[0]}=$row[1];
	$grid_coords_y{$row[0]}=$row[2];
    }
    close in;

    %patterns=();

    # Load all patterns in memory
    open (in,"$name$sw/$name$sw.patterns");
    while(<in>){
	chomp;
	@row=split(/\,/);
	$patterns{$row[0]}="";
	for my $loc (0..$locitot){
	    # Hash only the locations that are different from 
	    # the wild type "0" microsattelite length	    
	    # Shaka Zulu proposes hashing only the clone_id key
	    # and as value concatenate "locus,value,locus,value"
	    if($row[$loc+1]!=0){
		# Leading comma ",loc,val,loc,val"
		$patterns{$row[0]}=$patterns{$row[0]}.",".$loc.",".$row[$loc+1];
	    }
	}
    }
    close in;

    print STDERR "read done $sw\n";

    my $day = 182; 
    my $generation;
    my $firstline=1;
    # Go through the log file of grid state
    %grid=();
    open (i1,"$name$sw/$name$sw.log");
    open (out, ">$name$sw/$name$sw.divergence.log");
    while(<i1>){
	chomp;
	@row=split(/\,/);
	if($firstline){
            $generation=$row[3];
            $firstline=0;
        }
	if($row[2]!=$day){
	    
	    # compute divergence
	    #print STDERR "start comp diverg\n";
	    # Compute global divergence for $n cells
	    my $n = 300;
	    
	    # Draw cells
	    %hashxy = ();
	    %hashxy = draw_rand_cells($n,$grid_size_x);
	    %indexx = ();
	    %indexy = ();
	    my $i=0;
	    for my $x (keys %hashxy) {
		for my $y (keys %{$hashxy{$x}}) {
		    $indexx{$i}=$x;
		    $indexy{$i}=$y;
		    $i++;
		}
	    }
	    
	    $div = compute_divergence($n,\%grid,\%patterns,\%indexx,\%indexy);
	    $rounded = sprintf("%.6f", $div);
	    #print "$day,$rounded,$generation\n";
	    print out "$day,$rounded,$generation\n";
	    #print STDERR "diverg done\n";
	    # increase day
	    if($day==7280){
		$day=7300;
	    } else {
		$day+=182;
	    }
	    $generation=$row[3];
	    %grid=();
	    # grid X coordinate, Y coordinate, clone_id
	    #$grid{$row[0]}{$row[1]}=$row[2];
	    # Important - throw it as a double in hash
	    $grid{$grid_coords_x{$row[0]}*1.0}{$grid_coords_y{$row[0]}*1.0}=$row[1];
	} else {
	    # grid X coordinate, Y coordinate, clone_id
	    $grid{$grid_coords_x{$row[0]}*1.0}{$grid_coords_y{$row[0]}*1.0}=$row[1];
	    #print STDERR $grid_coords_x{$row[0]}." ".$grid_coords_y{$row[0]}." ".$row[1]."\n";
	}
    }
    close i1;    

    #print STDERR "loop done\n";

    # For day 7300
    # Compute global divergence for $n cells
    my $n = 300;
    
    # Draw cells
    %hashxy = ();
    %hashxy = draw_rand_cells($n,$grid_size_x);
    %indexx = ();
    %indexy = ();
    my $i=0;
    for my $x (keys %hashxy) {
	for my $y (keys %{$hashxy{$x}}) {
	    $indexx{$i}=$x;
	    $indexy{$i}=$y;
	    $i++;
	}
    }
    
    $div = compute_divergence($n,\%grid,\%patterns,\%indexx,\%indexy);
    #$div = compute_divergence($n);
    $rounded = sprintf("%.6f", $div);	    
    print out "$day,$rounded,$generation\n";
    close out;    
    print STDERR "computing divergence done $sw\n";
}

    
=pod
# Create allele dropouts
%new_patterns;

# Per panel 60% of 244 loci amplify
my %rem_hash = draw_rand_hash(147,244);
my %index_rem_hash, $ind=0;
for my $key (keys %rem_hash) {
    $index_rem_hash{$key} = $ind;
    $ind++;
}

for my $i (0..($n-1)){
    my $cid = $grid{$indexx{$i}}{$indexy{$i}};
    # Per individual cell, 60% of 147 loci amplify
    my %final_hash = draw_rand_hash(88,147);
    for my $loc (0..$locitot){
	$new_patterns{$cid}{$loc} = $patterns{$cid}{$loc};
	if(!$rem_hash{$loc}){
	    # Globally did not amplify
	    $new_patterns{$cid}{$loc} = -1;
	}
	if(!$final_hash{$index_rem_hash{$loc}}){
	    # Did not amplify in this cell
	    $new_patterns{$cid}{$loc} = -1;
	}
    }
}
=cut


sub compute_divergence{
    # Number of cells
    my ($n,$gridRef,$patternsRef,$indexxRef,$indexyRef) = @_;
    my %grid = %$gridRef;
    my %patterns = %$patternsRef;
    my %indexx = %$indexxRef;
    my %indexy = %$indexyRef;

    my $divergence = 0.0;
    my $dead_cells_count = 0;
    # Compare pattern of cell i to patterns of cell j
    # A total of n*(n-1)/2 comparisons 
    for my $i (0..($n-2)){
	my @col;
	# get the clone_id
	my $cid = $grid{$indexx{$i}}{$indexy{$i}};
#	print "CloneId ".$cid." to ".$cid2."\n";

	# Skip comparing dead cells
	if($cid==-1){
	    $dead_cells_count++;
	}
	# count if last cell is dead
	if($i==$n-2 && $grid{$indexx{$i+1}}{$indexy{$i+1}}==-1){
	    $dead_cells_count++;
	}

	if($cid!=-1){

	    for my $j (($i+1)..($n-1)){
		my $cid2 = $grid{$indexx{$j}}{$indexy{$j}};

		# Speedup: if cid2 is a dead cell, skip it.
		if($cid2!=-1){
		    # Significant speedup: If the two chosen cells are the same clone
		    # then divergence is 0 between them, thus skip looping over
		    # their microsattelite patterns
		    # Only when they are different clones, compare their pattern
		    if($cid!=$cid2){
			
			my $hamming = 0;
			#print "Cell ONE:\n";
			
			# Convert pattern
			my @cell_one = split(/\,/,$patterns{$cid});
			my @cell_two = split(/\,/,$patterns{$cid2});
			
			my %h1 = ();
			my %h2 = ();

			for my $crnch (1..$#cell_one){
			    if($crnch%2==1){
				# Hash locus, value
				$h1{$cell_one[$crnch]} = $cell_one[$crnch+1];
			    }
			}
			for my $crnch (1..$#cell_two){
			    if($crnch%2==1){
				# Hash locus, value
				$h2{$cell_two[$crnch]} = $cell_two[$crnch+1];
			    }
			}

			# Loop over cell one hashed pattern
			while (($locus, $value) = each(%h1)) {
			    #print "$cid $locus $value";
			    # If the locus is present in cell two
			    if(defined($h2{$locus})){
				$hamming += abs($h1{$locus}-$h2{$locus});			    
			    } else {
				$hamming += abs($h1{$locus});
			    }
			}
			#print "Cell TWO:\n";
			# Loop over cell two hashed pattern
			while (($locus, $value) = each(%h2)) {
			    #print "$cid2 $locus $value";
			    # If the locus is present in cell one, do nothing, if absent, compare to 100
			    if(!defined($h1{$locus})){			    
				$hamming += abs($h2{$locus});			 
			    }
			}
			# The rest of the loci are all 100s, so hamming distance is 0.
			# The distance between two cells is hamming/#of valid comparisons
			$divergence += $hamming*1.0/($locitot+1);
		    }
		}
	    }
	}
	#print "Comparing cell $i\n";
    }
    #print $n."\n";
    # At the end divide by the total comparisons between cells
    # subtract dead cells count from computation
    # What if all cells that we picked are dead? Then, set divergence to 0.
    if($dead_cells_count<$n-1){
	$divergence = $divergence/(($n-$dead_cells_count)*($n-$dead_cells_count-1)/2.0);
    } else {
	$divergence = 0.0;
    }
    return $divergence;
}

sub compute_divergence_filter{
    # Number of cells
    my ($n,$b1,$b2,$b3) = @_;
    #print "$b1 $b2 $b3\n";
    my $divergence = 0.0;
    # Compare pattern of cell i to patterns of cell j
    # A total of n*(n-1)/2 comparisons 
    my $pos=$b1;
    my $comp=0;
    for my $i (0..($n-2)){
	my @col;
	# get the clone_id
	my $cid = $grid{$indexx{$i}}{$indexy{$i}};
	# get the pattern
	for my $loc (0..$locitot){
	    $col[$loc]=$new_patterns{$cid}{$loc};
	}
	# bins

	if($i==$b1){
	    $pos=$b1+$b2;
	}
	if($i==$b1+$b2){
	    last;
	}
	for my $j ($pos..($n-1)){
	    $comp++;
	    #print "$comp ";
	    my $cid2 = $grid{$indexx{$j}}{$indexy{$j}};
	    #print "Grid ".$indexx{$i}." ".$indexy{$i}." to ".$indexx{$j}." ".$indexy{$j}."\n";
	    #print "CloneId ".$cid." to ".$cid2."\n";

	    my $state;
	    # zero total comparisons and hamming distance between two cells
	    my $total_comparisons = 0;
	    my $hamming = 0;
	    # get the pattern
	    for my $loc (0..$locitot){
		$state =$new_patterns{$cid2}{$loc};
		# if both states are present (neither states are -1)
		if($col[$loc]!=-1 && $state!=-1){
		    # count number of valid comparisons and hamming distance
		    $total_comparisons++;
		    $hamming += abs($col[$loc]-$state);
		    #print "States ".$col[$loc]." to ".$state."\n";
		}
	    }
	    # The distance between two cells is hamming/#of valid comparisons
	    $divergence += $hamming*1.0/$total_comparisons;
	}
	#print "Comparing cell $i\n";
    }
    #print $n."\n";
    # At the end divide by the total comparisons between cells
    $divergence = $divergence/($comp);
    return $divergence;
}






# Draw random x cells, but since we have hex grid
# x = 0, 0.5, 1, 1.5 ... 511, 511.5
# No two cells have the same coordinates
sub draw_rand_cells_limits{
    my ($n,$offs,$limxs,$limys) = @_;
    my %h=();
    my $size=0;
    while($size<$n){
        $x = int(rand(4));
	# either add 0.5 to both, or not
	$offset = 0.5*int(rand(2));
	$x+=$limxs+$offs+$offset;
	if($offset<0.2){
	    if($offs<0.2){
		$y = $limys+int(rand(10))*1.5;	    
	    } else {
		$y = $limys+int(rand(10))*1.5+0.75;	    
	    }
	} else {
	    if($offs<0.2){
		$y = $limys+int(rand(10))*1.5+0.75;	    
	    } else {
		$y = $limys+int(rand(10))*1.5;
	    }
	}
	# Ensure we pick out n DIFFERENT cells
	# do not pick out the same cell twice
        if(!$h{$x}{$y}){
            $h{$x}{$y}=1;
	    $size++;
        }
    }
    return %h;
}

=pod

# Get 100 x 3 biopsies from the segment
# to simulate real sampling
# Split the segment into three stripes and from each stripe get a biopsy
# For every 3 biopsies each 10 cells high and 40 cells wide
# take randomly from 5 to 30 cells
# and compute divergence
for my $biopsies (1..100){

    print "$biopsies,";
    # Choose coords of the biopsy
    # biopsy is 4mm by 1mm, 40 crypts by 10 crypts
    # 512 - 4
    my $x,$y,
    my $offset1;
    my $offset2;
    my $offset3;
    my $biopsy1xs, $biopsy1ys;
    my $biopsy2xs, $biopsy2ys;
    my $biopsy3xs, $biopsy3ys;
    # Bin1    
    #$x = int(rand(508));
    $x = int(rand($grid_size_x-4));
    $offset1 = 0.5*int(rand(2));
    $x += $offset1;
    $biopsy1xs = $x;
    if($offset1<0.2){
	# 85 - 10
	#$y = int(rand(75))*1.5;	    
	# 42-10
	# $grid_size_x/6 = 42
	$y = int(rand(32))*1.5;	    
	$biopsy1ys = $y;
    } else {
	$y = int(rand(32))*1.5+0.75;	    
	$biopsy1ys = $y;
    }	       	
    
    # Bin2
    $x = int(rand($grid_size_x-4));
    $offset2 = 0.5*int(rand(2));
    $x += $offset2;
    $biopsy2xs = $x;
    if($offset2<0.2){
	# 85 - 10
	$y = 42*1.5+int(rand(32))*1.5;	    
	$biopsy2ys = $y;
    } else {
	$y = 42*1.5+int(rand(32))*1.5+0.75;	    
	$biopsy2ys = $y;
    }	       	

    # Bin3
    $x = int(rand($grid_size_x-4));
    $offset3 = 0.5*int(rand(2));
    $x += $offset3;
    $biopsy3xs = $x;
    if($offset3<0.2){
	# 85 - 10
	$y = 84*1.5+int(rand(32))*1.5;	    
	$biopsy3ys = $y;
    } else {
	$y = 84*1.5+int(rand(32))*1.5+0.75;	    
	$biopsy3ys = $y;
    }	       	
    
 
    # for every biopsy, sample 3 to 30 cells
    for $cells (3..30){

	$n = $cells;
	# Choose how many cells to sample per bin, maximizing evennness
	my $bin1;
	my $bin2;
	my $bin3;
	if($cells%3==0){
	    $bin1=int($cells/3);
	    $bin2=int($cells/3);
	    $bin3=int($cells/3);	    
	} elsif($cells%3==1){
	    $bin1=int($cells/3);
	    $bin2=int($cells/3);
	    $bin3=int($cells/3)+1;    
	} else {
	    $bin1=int($cells/3);
	    $bin2=int($cells/3)+1;
	    $bin3=int($cells/3)+1;
	}
    
	my %hashbin1 = draw_rand_cells_limits($bin1,$offset1,$biopsy1xs,$biopsy1ys);

	my %hashbin2 = draw_rand_cells_limits($bin2,$offset2,$biopsy2xs,$biopsy2ys);

	my %hashbin3 = draw_rand_cells_limits($bin3,$offset3,$biopsy3xs,$biopsy3ys);

	%indexx=();
	%indexy=();
	my $i=0;
	
	# Combine hashbins
	#my %hashxy;
	for my $kx (keys %hashbin1){
	    for my $ky (keys %{$hashbin1{$kx}}){
		#$hashxy{$kx}{$ky}=1;
		$indexx{$i}=$kx;
		$indexy{$i}=$ky;
		$i++;	    
	    }
	}
	for my $kx (keys %hashbin2){
	    for my $ky (keys %{$hashbin2{$kx}}){
		#$hashxy{$kx}{$ky}=1;
		$indexx{$i}=$kx;
		$indexy{$i}=$ky;
		$i++;	    
	    }
	}
	for my $kx (keys %hashbin3){
	    for my $ky (keys %{$hashbin3{$kx}}){
		#$hashxy{$kx}{$ky}=1;
		$indexx{$i}=$kx;
		$indexy{$i}=$ky;
		$i++;
	    }
	}

	# Draw cells within boundaries of biopsy
	#my %hashxy = draw_rand_cells($n,512);
	#for my $x (keys %hashxy) {
	#    for my $y (keys %{$hashxy{$x}}) {
		#$indexx{$i}=$x;
		#$indexy{$i}=$y;
		#$i++;
	    #}
	#}
    
	# Create allele dropouts
	%new_patterns=();
    
	# Per panel 60% of 244 loci amplify
	my %rem_hash=draw_rand_hash(147,244);
	my %index_rem_hash=();
	my $ind=0;
	for my $key (keys %rem_hash) {
	    $index_rem_hash{$key} = $ind;
	    $ind++;
	}
    
	for my $i (0..($n-1)){
	    my $cid = $grid{$indexx{$i}}{$indexy{$i}};
	    # Per individual cell, 60% of 147 loci amplify
	    my %final_hash = draw_rand_hash(88,147);
	    for my $loc (0..$locitot){
		$new_patterns{$cid}{$loc} = $patterns{$cid}{$loc};
		if(!$rem_hash{$loc}){
		    # Globally did not amplify
		    $new_patterns{$cid}{$loc} = -1;
		}
		if(!$final_hash{$index_rem_hash{$loc}}){
		    # Did not amplify in this cell
		    $new_patterns{$cid}{$loc} = -1;
		}
	    }
	}
    
	#print "Computing divergence for $cells cells...\n";
	#$div = compute_divergence($n);
	$div = compute_divergence_filter($n,$bin1,$bin2,$bin3);
	if($cells==30){
	    $rounded = sprintf("%.6f", $div);
	    print "$rounded\n";	    
	} else {
	    $rounded = sprintf("%.6f", $div);
	    print "$rounded,";
	}
    }    
}
















=cut

